# Import Packages
pacman::p_load(sf, tmap, tidyverse, sfdep, readxl, Kendall, plotly, plyr)EDA
# Import MPSZ
mpsz <- st_read(dsn = "data/geospatial", layer = "MPSZ-2019")Reading layer `MPSZ-2019' from data source
`C:\wamp64\www\IS415-Project\data\geospatial' using driver `ESRI Shapefile'
Simple feature collection with 332 features and 6 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 103.6057 ymin: 1.158699 xmax: 104.0885 ymax: 1.470775
Geodetic CRS: WGS 84
mpsz <- st_transform(mpsz, 3414)
st_crs(mpsz)Coordinate Reference System:
User input: EPSG:3414
wkt:
PROJCRS["SVY21 / Singapore TM",
BASEGEOGCRS["SVY21",
DATUM["SVY21",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4757]],
CONVERSION["Singapore Transverse Mercator",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",1.36666666666667,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",103.833333333333,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",1,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",28001.642,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",38744.572,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["northing (N)",north,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["easting (E)",east,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Cadastre, engineering survey, topographic mapping."],
AREA["Singapore - onshore and offshore."],
BBOX[1.13,103.59,1.47,104.07]],
ID["EPSG",3414]]
write_rds(mpsz, "data/model/mpsz.rds")# Import MRT Stations
mrt <- read.csv("data/geospatial/mrtsg.csv")
mrt_sf <- st_as_sf(mrt,
coords = c("Longitude",
"Latitude"),
crs = 4326) %>%
st_transform(crs = 3414)write_rds(mrt_sf, "data/model/mrt_sf.rds")# Tourism
tourism_sf <- st_read(dsn = "data/geospatial", layer = "TOURISM")Reading layer `TOURISM' from data source
`C:\wamp64\www\IS415-Project\data\geospatial' using driver `ESRI Shapefile'
Simple feature collection with 107 features and 19 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: -1.797693e+308 ymin: -1.797693e+308 xmax: 43659.54 ymax: 47596.73
Projected CRS: SVY21
# Assign EPSG Code
tourism_sf <- st_transform(tourism_sf, 3414)write_rds(tourism_sf, "data/model/tourism_sf.rds")# Shopping Malls
shopping <- read.csv("data/geospatial/mall_coordinates.csv")
shopping <- shopping %>%
select(name, latitude, longitude)
glimpse(shopping)Rows: 199
Columns: 3
$ name <chr> "100 AM", "i12 Katong", "313@SOMERSET", "321 CLEMENTI", "600…
$ latitude <dbl> 1.274588, 1.305087, 1.301385, 1.312025, 1.334042, 1.437131, …
$ longitude <dbl> 103.8435, 103.9051, 103.8377, 103.7650, 103.8510, 103.7953, …
shopping_sf <- st_as_sf(shopping,
coords = c("longitude",
"latitude"),
crs = 4326) %>%
st_transform(crs = 3414)write_rds(shopping_sf, "data/model/shopping_sf.rds")# Hexagon for Map
hexagons <- st_read(dsn = "data/geospatial", layer = "hexagons")Reading layer `hexagons' from data source
`C:\wamp64\www\IS415-Project\data\geospatial' using driver `ESRI Shapefile'
Simple feature collection with 3125 features and 6 fields
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: 2667.538 ymin: 21506.33 xmax: 50010.26 ymax: 50256.33
Projected CRS: SVY21 / Singapore TM
hexagons_sf <- st_transform(hexagons, 3414)write_rds(hexagons_sf, "data/model/hexagons_sf.rds")# HDB Population by Geographical Location
hdb <- read.csv("data/aspatial/hdb-resident-population-by-geographical-distribution.csv")# Filter out 2018 data
hdb <- filter(hdb, shs_year == "2018")write_rds(hdb, "data/model/hdb.rds")tmap_mode('view')
tmap_options(check.and.fix = TRUE) +
tm_shape(mpsz) +
tm_polygons() +
tm_shape(mrt_sf) +
tm_dots(alph=0.5, size=0.01)+
tm_view(set.zoom.limits = c(11,14))mpsz_hdb <- left_join(mpsz, hdb, by = c("PLN_AREA_N" = "town_estate"))# Population Density
tmap_mode('view')
tm_shape(mpsz_hdb) +
tm_fill("number")# MRT vs Population Density
tmap_mode('view')
tm_shape(mpsz_hdb) +
tm_fill("number") +
tm_shape(mrt_sf) +
tm_dots(alph=0.5, size=0.01)+
tm_view(set.zoom.limits = c(11,14))# Population Data 2022
popdata <- read_csv("data/aspatial/respopagesexfa2022.csv")popdata2022 <- popdata %>%
spread(AG, Pop) %>%
mutate(YOUNG = `0_to_4`+`5_to_9`+`10_to_14`+
`15_to_19`+`20_to_24`) %>%
mutate(`ECONOMY ACTIVE` = rowSums(.[9:13])+
rowSums(.[15:17]))%>%
mutate(`AGED`=rowSums(.[18:22])) %>%
mutate(`TOTAL`=rowSums(.[5:22])) %>%
mutate(`DEPENDENCY` = (`YOUNG` + `AGED`)
/`ECONOMY ACTIVE`) %>%
mutate_at(.vars = vars(PA, SZ), toupper) %>%
select(`PA`, `SZ`, `YOUNG`, `ECONOMY ACTIVE`, `AGED`,
`TOTAL`, `DEPENDENCY`) %>%
filter(`ECONOMY ACTIVE` > 0)popdata2022 <- popdata2022 %>%
mutate_at(.vars = vars(PA, SZ),
.funs = funs(toupper)) %>%
filter(`ECONOMY ACTIVE` > 0)write_rds(popdata2022, "data/model/popdata2022.rds")mpsz_popdata2022 <- left_join(mpsz, popdata2022, by = c("SUBZONE_N" = "SZ"))write_rds(mpsz_popdata2022, "data/model/mpsz_popdata2022.rds")# MRT vs Population Density
tmap_mode('view')
tm_shape(mpsz_popdata2022) +
tm_fill("TOTAL") +
tm_borders(lwd = 0.1, alpha = 1) +
tm_shape(mrt_sf) +
tm_dots(alph=0.5, size=0.01)+
tm_view(set.zoom.limits = c(11,14))